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Abstract 

Finite element analysis has been used successfully to estimate the effective properties of many types of composites. 
The prediction of effective elastic moduli of polymer-bonded explosives provides a new challenge. These particulate 
composites contain extremely high volume fractions of explosive particles (> 0.90). At room temperature and higher, 
the Young's modulus of the particles can be 20,000 times that of the binder. Under these conditions, rigorous bounds 
and analytical approximations for effective elastic properties predict values that are orders of magnitude different 
from the experimental values. In this work, an approach is presented that can be used to predict three-dimensional 
effective elastic moduli from two-dimensional finite element simulations. The approach is validated by comparison 
with differential effective medium estimates and three-dimensional finite element simulations. The two-dimensional 
finite element approach has been used to determine the properties of models of polymer bonded explosives and PBX 
9501 in particular, containing high volume fractions of circular and square particles with high modulus contrasts. 
Results show that estimates of effective elastic properties from two-dimensional finite element calculations are close 
to the values predicted by the differential effective medium approach for a large range of volume fractions and 
modulus contrasts. Two- and three-dimensional finite element estimates for volume fractions from 0.70 to 0.90 and 
found not to differ considerably. Simulations of models of polymer bonded explosives and PBX 9501 show that the 
microstructure, the amount of discretization, and the type of element used play a considerable role in determining the 
value predicted by finite element simulations. The effective elastic moduli of PBX 9501 predicted by finite element 
calculations can vary from 200 MPa to 10,000 MPa depending on the microstructure and level of discretization used. 
The results also suggest that if a microstructure can be found that accurately predicts the elastic properties of PBX 
9501 at a certain temperature and strain rate, then the same microstructure can be used to predict elastic properties at 
other temperatures and strain rates. 

Keywords : particle-reinforced composites, microstructure, modeling, elastic properties 

1 Introduction 

Experimental determination of the mechanical properties of polymer-bonded explosives (PBXs) is hazardous and 
highly expensive. An alternative to experimentation is the determination of material properties using micromechanics 
methods. However, rigorous bounds and analytical approaches predict values of elastic modulus that are orders of 
magnitude different from the experimental values. Hence, detailed numerical simulations are required for the 
prediction of effective properties of PBXs. The goal of this work is to investigate if the elastic properties of PBXs can 
be determined using finite element analyses of PBX microstructures. 

*E-mail: banerjee@eng.utah.edu Fax: (801)-585-00396 
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PBXs are particulate composites containing explosive particles and a rubbery binder The particle size distributions in 
these materials and the mechanical properties of the constituents are required prior to finite element simulations of 
the microstructure of these materials. The size distribution of the explosive particles in PBXs can be obtained 
experimentally. The mechanical properties of the particles can be determined from molecular dynamics simulations. 
Mechanical testing of the rubbery binder does not involve any risk of explosion. Hence, mechanical properties of the 
binder can be determined from laboratory tests. 

The major challenges involved in finite element modeling of PBXs are the high volume fraction of particles in PBXs 
(/p > 0.90) and the high modulus contrast between particles and binder (Ep/Eb = 10,000 to 20,000) at and above 
room temperature and at low strain rates. Because of the high volume fractions, microstructures are difficult to 
generate and even more difficult to digitize from micrographs. In addition, the high modulus contrast leads to large 
stress concentrations at the interface of particles and the binder. These stress concentrations cannot be simulated 
adequately without a high degree of mesh discretization. Hence, three-dimensional finite element models of these 
materials require considerable computational power. 

In this work, a procedure is outlined whereby two-dimensional finite element simulations can be used to estimate the 
three-dimensional elastic moduli. In order to verify the accuracy of the approach, finite element estimates on unit 
cells containing circular particles are compared with differential effective medium estimates Uji2,| of effective elastic 
moduli for a wide range of volume fractions and modulus contrasts. Two-dimensional finite element estimates for 
selected microstructures are then compared to three-dimensional estimates for particle volume fractions of 0.70, 0.75 
and 0.80. 

Next, an appropriate model of PBX 9501 is determined from simulations of manually generated and randomly 
generated microstructures based on the particle size distribution of PBX 9501. Finite element calculations are then 
performed on this model to estimate the elastic moduli of PBX 9501 for various temperatures and strain rates. 
Finally, the estimated moduli are compared with experimental data to gauge the applicability of the approach to 
polymer bonded explosives. 



2 Finite element approach 



Numerical determination of the effective properties of particulate composites involves the calculation of the stress 
and strain fields for a representative volume element (RVE) that simulates the microstructure of the composite. These 
stresses and strains are averaged over the volume (V) of the RVE. The effective elastic stiffness tensor C|^; of the 
composite is calculated from the tensor relation 
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where cTy are the stresses and are the strains. 

The stress-strain relation for the RVE can be written in Voigt notation as 
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where {a)y is the volume average of the quantity a. Equation Q can also be inverted and written in terms of the 
directional Young's moduli, Poisson's ratios and shear moduli as 
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where the directional Young's moduli Ef[, E^, and differ by a small amount from the isotropic Young's 
modulus i?eff. Similarly, the mean of the Poisson's ratios i^^f , v?^, vf^, , j^lf , and i^§f is assumed to be close to the 
isotropic Poisson's ratio v^a and the mean of the shear moduli G^f , GI3, and is assumed to be close to the 
isotropic shear modulus Geff. The bulk mechanical response of a random particulate composite is isotropic. However, 
the size of a RVE of such a composite that can be simulated numerically is necessarily much smaller than the bulk 
and may not be isotropic. In this work, it is assumed that the deviation from isotropy of an RVE is small. 
The determination of effective isotropic elastic modulus of particulate composites using three-dimensional finite 
element analysis therefore requires the determination of volume averaged stresses and strains under appropriate 
boundary conditions. Once these stresses and strains are known for an RVE, equation (|3} can be used to determine 
the effective elastic moduli of the RVE. An investigation into effects of boundary conditions revealed a uniaxial state 
of stress develops in a RVE when a uniform boundary displacement is applied in one direction and periodicity of 
boundary displacements is maintained. This means that the effective Young's modulus and Poisson's ratio in a 
particular direction can be directly computed from equation Q by applying a uniform displacement in that direction. 
The effective isotropic Young's modulus and Poisson's ratio can then be calculated by averaging the values obtained 
in the three orthogonal directions from three simulations. The effective shear modulus can be calculated directly from 
the effective isotropic Young's modulus and Poisson's ratio. 

The estimation of effective elastic properties from three-dimensional simulations is straightforward in principle. 
However, not only is it difficult to generate RVEs containing high volume fractions of random spherical particles, 
meshing such three-dimensional RVEs is nontrivial. Even if RVEs can be generated and meshed, solution of the finite 
element system of equations can be prohibitive in terms of computational cost. A two-dimensional approach is 
described below that requires relatively little computational expense. This approach can be used to arrive at 
reasonable estimates of three-dimensional effective properties of particulate composites. 
In the two-dimensional simulations of particulate composites performed in this study, a cross-section of the 
composite is chosen that contains the same area fraction of particles as the volume fraction of particles in the 
composite. Finite element analyses are then performed on this square RVE assuming a state of plane strain. The 
stress-strain equation (|3} can be written in planar or two-dimensional form as 



(eii)v 




(e22)v 




2(ei2)v 





, ,eff / Rieff 
-f^l2/-'i22 














(c^ii)v 







(C^22)v 






>12)v. 



(4) 



where Eff, ii^lf , I'fl, and z^jf are two-dimensional Young's moduli and Poisson's ratios and different from their 
three-dimensional counterparts. A brief explanation of two-dimensional elastic moduli is given in AppendixlAl 
Uniform displacements are applied to the '1' and '2' directions of the RVE and periodic displacements are enforced 
along the remainder of the boundary to arrive at uniaxial average stresses and the corresponding average strains in the 
RVE. The effective moduli are calculated using equation (0} and the directional moduli are averaged to arrive at a 
mean effective Young's modulus (E^^) and a mean effective Poisson's ratio (I'g^) for the RVE. These moduli can 
also be thought of as effective transverse moduli of a composite containing unidirectional cylinders. 
It is assumed that the transverse moduli are equal to the two-dimensional apparent moduli of an isotropic, 
homogeneous material subjected to plane strain. Since equation is a two-dimensional stress-strain law, the 
two-dimensional moduli have to be converted to three-dimensional moduh using the relations Q 



E,s = i?,2t?[l-(i^eff)'], 



(5) 



where v^ff is the effective Poisson's ratio in three dimensions and E^ff is the three-dimensional effective Young's 
modulus. 
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3 Validation of approach 



3.1 Two-dimensional FEM vs. differential effective medium estimates 

The differential effective medium (DEM) approach (fllJl has been found to provide excellent estimates of elastic 
properties of particulate composites (3- To determine if the two-dimensional finite element approach provides 
reasonable estimates of three-dimensional effective elastic properties, finite element analyses were performed on the 
two-dimensional RVEs containing 10% to 92% by volume of circular particles that are shown in Figure^ The 
particles were assigned a Young's modulus of 100,000 MPa and a Poisson's ratio of 0.2. The Young's modulus of the 
binder was varied from 1 MPa to 10,000 MPa in multiples of 10 and the binder Poisson's ratio was chosen to be 0.49. 
The effective properties predicted by the finite element approach were then compared with DEM estimates. 
In the RVEs shown in Figure[n circular particles were placed sequentially at random locations in each RVE such that 
no two particles were in contact. The particle size distribution in each RVE was based on that of the dry blend of 
PBX 9501 |4|. The planar moduli obtained from the finite element calculations were converted into 
three-dimensional moduli using equations|5] Figure|2ja) shows effective three-dimensional Young's moduli of the 
RVEs from finite element (FEM) and differential effective medium (DEM) calculations. Figure|2jb) shows the 
corresponding three-dimensional Poisson's ratios. The solid lines in the figures are lines of constant modulus contrast 
between particles and binder and represent the effective elastic property calculated using DEM. The finite element 
results are shown as squares. 

The FEM and the DEM predictions are seen to agree remarkably well for all modulus contrasts and for all the 
simulated volume fractions of particles, implying that the two-dimensional FEM approach described in the previous 
section gives reasonable estimates of effective elastic properties over a considerable range of volume fractions and 
modulus contrasts. 

3.2 Two-dimensional vs. three-dimensional FEM estimates 

To obtain a more direct estimate of the effectiveness of the two-dimensional FEM approach, a second set of finite 
element simulations was performed on two- and three-dimensional RVEs containing particle volume fractions of 
0.70, 0.75, and 0.80. The two-dimensional RVEs are shown in Figure|3ja) and the three-dimensional RVEs are shown 
in Figure|3b). The Young's modulus used for the particles was 15,300 MPa and the Poisson's ratio was 0.32. The 
Young's modulus of the binder was 0.7 MPa and the Poisson's ratio was 0.49. These values correspond to the moduli 
of HMX and the binder of PBX 9501 at room temperature and low strain rate |5 6 1. 

For the two-dimensional simulations, each RVE was discretized using six-noded triangles and approximately 60,000 
nodes. Displacement boundary conditions were applied and the average stresses and strain were used to calculate 
two-dimensional effective Young's moduli and Poisson's ratios. The two-dimensional moduli were then converted 
into three-dimensional moduli using equations|5] The three-dimensional RVEs were discretized using 10-noded 
tetrahedral elements with a total of around 100,000 nodes. A uniform displacement was applied perpendicular to one 
of the faces of the RVE while the remaining faces were constrained to displace in a periodic manner. The 
three-dimensional Young's modulus and Poisson's ratio were directly calculated from the average stresses and strains 
for these three-dimensional models. 

FigureEJa) shows the predicted effective Young's moduli of the two- and three-dimensional RVEs shown in Figure|3 
Though the three-dimensional models predict higher values of Young's modulus than the two-dimensional models, 
the difference between the two- and three-dimensional estimates of Young's modulus is small relative to the modulus 
contrast between the particles and the binder. Figure01b) shows the effective Poisson's ratio for the two- and 
three-dimensional finite element models. The three-dimensional models predict lower values of Poisson's ratio than 
the two-dimensional models. Both two- and three-dimensional models predict a sharp decrease in Poisson's ratio 
with increasing particle volume fraction. These results show that the two-dimensional FEM approach can be used to 
obtain acceptable estimates of effective elastic properties of random particulate composites and hence as a substitute 
for detailed three-dimensional calculations. 
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4 Modeling polymer bonded explosives 



A micrograph of the polymer bonded explosive PBX 9501 Q is shown in Figure|5l The HMX particles are shown to 
be irregularly shaped and present in a large number of sizes. Two length scales can be identified from the micrograph. 
The first is the scale of the larger particles that occupy most of the volume. The second scale is that of the particles 
filling the interstitial spaces between the larger particles. Because of these two different length scales, it is extremely 
difficult to use a digital image at a single scale to generate two-dimensional microstructures of PBX 9501 or other 
PBXs. 

Most micromechanical calculations for PBXs have been carried out using subgrid models that use simplified models 
of the microstructure (for example, spherical grains coated with binder or spherical voids in an effective PBX 
material 1 8 1). Closed form solutions from these simple models have been used to provide properties for macroscopic 
simulations. More detailed calculations have used microstructures containing ordered arrays of circles or polygons in 
two or three dimensions to model PBXs |9 10 1. These models do not reflect the microstructure of PBXs and hence 
have limited use for predicting thermoelastic properties. Better two-dimensional approximations of the 
microstructure have been constructed from digital images of the material and used by Benson and Conley 11 II to 
study some aspects of the micromechanics of PBXs. However, such microstructures are extremely difficult to 
generate and require state-of-the-art image processing techniques and excellent images to accurately capture details 
of the geometry of PBXs. More recently Baer 1 12 1 has used a combination of Monte Carlo and molecular dynamics 
techniques to generate three-dimensional microstructures containing spheres and oriented cubes that appear to 
represent PBX microstructures well. However, the generation of a single realization of these microstructures is very 
time consuming and often leads to a maximum packing fraction of about 70%. Periodicity is also extremely difficult 
to maintain in the RVEs generated by this method. 

One of the goals of this work was to explore a fast, yet accurate, means of predicting the effective elastic properties of 
PBXs. With this goal in mind, accurate representation of the microstructure of PBXs was foregone. Instead, the 
explosive crystals were modeled using circles or squares. 

Both manually and automatically generated two-dimensional RVEs of PBX 9501 are discussed in this section. The 
effective elastic moduli of these RVEs from finite element calculations are then compared with experimentally 
determined moduli of PBX 9501. As in the previous section, the particles in the RVEs represent HMX and have a 
Young's modulus of 15,300 MPa and a Poisson's ratio of 0.32 |5 1 and the binder has a Young's modulus of 0.7 MPa 
and a Poisson's ratio of 0.49 ||6l. The effective properties computed from these component properties are compared 
with flie Young's modulus of 1,013 MPa and Poisson's ratio of 0.35 1^61 of PBX 9501. 



4.1 Manually generated microstructures 

Figure |6l shows six two-dimensional microstructures representing PBX 9501. Each RVE is filled with particles of 
various sizes since PBXs are typically a mixture of coarse and fine grains with the finer grains forming a filler 
between coarser grains. The volume fractions of particles in each of these models is 90±0.5%. All the models 
possess square symmetry. The RVEs have also been designed so that particle-particle contact is avoided. Finite 
element simulations were performed on these RVEs using six-noded triangular elements such that the circular 
particles were discretized accurately. The approach discussed in Section|2lwas used to determine the effective 
Young's modulus and Poisson's ratios of these RVEs. 

Table^shows the effective Young's modulus and Poisson's ratios of the six RVEs. The values of the effective 
Young's modulus range from 42 MPa to 192 MPa, with a mean of 132 MPa and a standard deviation of 54 MPa. In 
comparison, the Young's modulus of PBX 9501 is 1,013 MPa; around 500% to 800% of that predicted using the six 
RVEs. Among the six RVEs, the Young's modulus appears to be higher for the models with lower amounts of binder 
along the edges of the RVE. Models 1 through 3 have a single large particle and many smaller particles and show 
approximately the same effective behavior Model 4, with a smaller ratio between the radius of the largest and the 
smallest particles, has a lower effective Young's modulus than the mean value. 

The effective Poisson's ratios of the six RVEs range from 0.25 to 0.34, with a mean of 0.33. The mean Poisson's ratio 
of the RVEs is not much different from the measured value for PBX 9501. For models 1, 2, and 3, the effective 
Poisson's ratio is within 10% of that of PBX 9501. However, the Poisson's ratio of model 4 is 25% higher and those 
of models 5 and 6 are 20% and 28% lower than the Poisson's ratio of PBX 9501, respectively. In comparison, the 
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two-dimensional Poisson's ratios of models 4, 5, and 6 are 0.79, 0.39, and 0.33, respectively, whereas that of PBX 
9501 is 0.54. Since these two-dimensional values reflect the actual strains obtained from the two-dimensional finite 
element calculations, it is clear that when compared to PBX 9501, model 4 is too compliant while models 5 and 6 are 
too stiff. These observations are also attributable to the amount of binder along the edge of a RVE. 
On average, the FEM estimate of Young's modulus is around 13% of the experimentally determined value for PBX 
9501. The reason for this difference could be that preferential stress paths (stress-bridging) in the microstructure of 
PBX 9501 lead to increased stiffness. Results shown in Figurel^indicate that the difference between estimates from 
two- and three-dimensional models is not large. Hence, it is unlikely that the considerably lower stiffness of the six 
model RVEs compared to PBX 9501 is due to the use of two-dimensional models. 

It is possible that the Young's modulus is underestimated by the six RVEs because the volume fraction of particles in 
each of these microstructures is lower than 92%. Figure0a) shows a model containing 92% particles by volume that 
is a modified version of model 6 in Figure|6l The effective Young's modulus of the RVE containing 92% particles is 
218 MPa and the effective Poisson's ratio is 0.28. In comparison, model 6 contains around 90% particles by volume 
and has an estimated Young's modulus of 192 MPa and a Poisson's ratio of 0.25. Though there is some increase in 
the Young's modulus due to increase in particle volume fraction, the value is still around 20% that of PBX 9501. 
Hence, the volume fraction cannot be the only factor leading to underestimation of the Young's modulus by the 
manually generated microstructures. 

If a 256x256 square grid is overlaid on top of the model shown in Figure0a) and cells in the grid are assigned HMX 
properties if they contain more than 50% particle and binder properties otherwise (referred to as the "square grid 
overlay method"), the new microstructure shown in Figure0d) is obtained. The RVE in FigureQd) was modeled 
using 256 X 256 four-noded square finite elements. The modified model predicts a Young's modulus of 800 MPa and 
a Poisson's ratio of 0.14. The modulus is now closer to that of PBX 9501 but the Poisson's ratio is much lower. 
The stiffer response of the model shown in Figure0d) is partially due to increased stiffness of four-noded finite 
elements. To eliminate the possibility that the stiffer response of the modified model is due to element locking caused 
by the nearly incompressible binder, the binder was also modeled as a two-parameter Mooney-Rivlin rubber instead 
of a linear elastic material. An incremental analysis using the Mooney-Rivlin material model yielded values of 
Young's modulus and Poisson's ratio that were within 5% of those for the linear elastic models. This result implies 
that the difference between the effective properties of the RVEs in FiguresQc) and^Id) is primarily due to the slight 
difference in microstructure. The same behavior has also been observed for all the models shown in Figure|6l 
The model in FigureQd) has a Young's modulus that is four times higher than that of the model in Figure0c). This 
result suggests that the square grid overlay method can lead to the automatic incorporation of preferential stress paths 
into RVEs containing high volume fractions of particles. In addition, it is easier to generate a square mesh and use the 
square grid overlay method than to model circular particles using triangular elements. Hence, the square grid overlay 
method is used to discretize the randomly generated microstructures based on the actual particle size distribution of 
PBX 9501 discussed in the next section. 

The manually generated microstructures were based on a qualitative assessment of the particle size distribution of 
PBX 9501 without regard for the actual size distribution. In addition, particles in the RVEs were not allowed to come 
into direct contact leading to a gross underestimation of the effective Young's modulus. The first of these issues is 
addressed in the following section where the actual particle size distribution of PBX 9501 is used to generate model 
microstructures. The particle contact problem is addressed by approximating the model microstructures using the 
square grid overlay method. 

4.2 Randomly generated microstructures 

The preferred method for generating close packed microstructures from a set of particles is to use Monte Carlo based 
molecular dynamics techniques |T3l or Newtonian motion based techniques 1 14|. In both methods, a distribution of 
particles is allocated to the grid points of a rectangular lattice using a random placement method. Molecular 
dynamics simulations or Newtonian dynamics calculations are then carried out on the system of particles to reach the 
packing fraction that corresponds to equilibrium. A weighted Voronoi tessellation 1 15 1 is then carried out on the 
particles with the weights determined by the sizes of the particles. The particles are next moved towards the center of 
the packing volume while maintaining that they remain inside their respective Voronoi cells. The process is repeated 
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until all the particles are as tightly packed as possible. Periodicity of the particles at the boundaries is maintained by 
specifying extra particles at the boundaries that move in and out of the volume. This process, with some 
modifications, has been the only efficient method of generating close packed systems of particles in three dimensions. 
However, it is difficult to get packing fractions of more than 70-75% when using spherical particles. It is also quite 
time consuming to generate a tight packing. 

In two dimensions, a faster method can be used for generating particle distributions, that is random sequential 
packing. The largest particles are placed randomly in the volume followed by progressively smaller particles. If there 
is any overlap between a new particle and the existing set, the new particle is moved to a new position. If a particle 
cannot be placed in the volume after a certain number of iterations, the next lower sized particle is chosen and the 
process is continued until the required volume fraction is achieved. Though this method does not preserve the particle 
size distributions as accurately as the Monte Carlo based molecular dynamics methods, it is much faster and can be 
used to generate high packing fractions in two dimensions without particle locking. In three dimensions, this method 
is highly inefficient and packing above 60% is extremely time consuming to achieve. 



4.2.1 Circular particles - PBX 9501 dry blend 

HMX particle size distributions of the dry blend of PBX 9501 have been listed by Wetzel JS) and Rae et al. 1 16 1. The 
coarse and fine particles are blended in a ratio of 1:3 by weight and compacted to form PBX 9501. Four 
microstructures based on the particle size distribution of the dry blend are shown in Figure|8l The number of particles 
used for the four microstructures are 100, 200, 300, and 400. The RVE widths are 0.65 mm, 0.94 mm, 1.13 mm and 
1.325 mm, respectively. The particles occupy a volume fraction of about 85-86%. The remaining volume is assumed 
to be occupied by a mixture of binder and fine particles of HMX that are well separated in size from the smallest 
particles shown in the RVEs. This 'dirty' binder contains around 36% particles by volume. 
Particles in the RVEs were assigned HMX properties and the elastic properties of the dirty binder were calculated 
using the differential effective medium approximation The Young's modulus of the dirty binder used in the 
calculations was 2.0816 MPa and the Poisson's ratio was 0.4813. Two sizes of square grids were used for overlaying 
on the RVEs: 256x256 and 350x350. The square grid overlay method was applied to assign materials to grid cells. 
Each cell in the grid was modeled using a four-noded finite element. Table|2]shows the effective Young's modulus 
and Poisson's ratio predicted for the four models of the dry blend of PBX 9501. 

The FEM calculations overestimate the effective Young's moduh of PBX 9501 by 200% to 400% if the 256x256 
grid is used. If the 350 x 350 grid is used the effective Young's moduH vary from 100% to 300% of that of PBX 9501 . 
The effective Poisson's ratio is underestimated considerably in both cases. The stiffer response of these RVEs is 
partly due to the creation of continuous stress paths across a RVE when a regular grid is overlaid on the RVE and 
partly because four-noded elements are inherently stiffer [171. It can also be seen that the effective Young's modulus 
increases as the number of particles in the RVE increases. This is because the same grid has been used for the four 
microstructures, leading to poorer resolution of the geometry with increasing particle density per grid cell. 



4.2.2 Circular particles - pressed PBX 9501 

Figure|9]shows four RVEs based on the particle size distribution of pressed PBX 9501 fT^. The pressing process leads 
to particle breakage and hence a larger volume fraction of smaller sized particles. Fewer larger sized particles remain 
as is reflected in the generated microstructures containing 100, 200, 500, and 1000 particles. In this case, the RVE 
widths are 0.36 mm, 0.42 mm, 0.535 mm, and 0.68 mm, respectively. Thus, the 1000 particle RVE for the pressed 
piece has dimensions similar to the 100 particle RVE for the dry blend. The size of the RVE that can be adequately 
discretized is therefore smaller for the pressed piece. Each RVE was discretized into both 256x256 and 350x350 
four-noded elements. Each element was assigned material properties according to the square grid overlay method. 
The particles in the RVEs occupy volume fractions from 0.86 to 0.89 and the target volume fraction of 0.92 is not 
attained in any of the RVEs. A dirty binder, whose properties were calculated using the differential effective medium 
approach, was used in the effective stiffness calculations. The effective moduli of the dirty binder for the four models 
are shown in Table|3] The effective Young's modulus and Poisson's ratio of the four RVEs from FEM calculations are 
shown in Table|4] The 256 x256 element models overestimate the Young's modulus of PBX 9501 by factors 
increasing from 3 to 6 with increasing RVE size. The estimates from the 350 x 350 element models are lower but still 
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2 to 5 times higher than the Young's modulus of PBX 9501 . The estimated Poisson's ratios are quite low compared to 
thatofPBX 9501. 

The four models based on pressed PBX 9501 predict Young's moduli that are 1.5 to 2 times higher than values 
predicted by models based on the dry blend of PBX 9501. For the 100 and 200 particle pressed PBX 9501 models, 
the single large particle contributes considerably to the stiffer response. For the 500 and 1000 particle pressed PBX 
9501 models, errors in the discretization of particle boundaries lead to additional stress bridging paths and hence a 
stiffer response is obtained. 

4.2.3 Square particles - pressed PBX 9501 

With an increase in particle volume fraction, there is an increase in the number of particle-particle contacts in 
particulate composites. If the particles in a RVE of such a material are circular, triangular elements cannot be used to 
discretize the RVE because of poor element shapes in regions of contact. Additionally, a rectangular grid cannot 
represent the geometry of circular particles accurately. An alternative is to use square particles that are aligned with a 
rectangular grid to represent the microstructure. 

The particles shown in the three microstructures in Figure[^are based on the size distribution of pressed PBX 9501. 
These distributions have been generated by placing square particles in a random sequential manner in a square grid. 
The smallest particles occupy a single subcell of the grid. Larger particles are chosen from the particle size 
distribution so that they fit into an integer multiple of the grid size. 

In the three models shown in FigurefTO] the particle size distribution for the pressed piece is truncated so that the 
smallest (< 30 microns) particle sizes in the distribution are not used in order that the grid size does not become too 
large. The RVEs are filled with particles to volume fractions of about 86-87%. The remaining volume is assumed to 
be occupied by a dirty binder The number of particles in the first model is 700 and the RVE width is around 3.6 mm. 
The second model contains 2,800 particles and the RVE width is about 5.3 mm, while the third model contains 
1 1,600 particles and has a width of 9 mm. Note that these RVEs are considerably larger than those used for the 
circular particles. Smaller RVEs are not used because of the difficulties associated with fitting particles into integer 
multiples of subcell widths. The properties of the dirty binder for the three models are shown in Table|5] 
Finite element calculations were performed on the three models using regular grids of 256x256 square elements. 
Table|6lshows the effective Young's modulus and Poisson's ratio obtained from the three RVEs. The FEM estimates 
of Young's modulus are 9 times higher than the Young's modulus of PBX 9501 . The higher stiffness obtained from 
the FEM calculations is partly because the large number of contacts between particles lead essentially to a continuous 
particle phase containing pockets of binder. 

Micrographs of PBX 9501 show that there is considerable damage in the particles, binder and the particle-binder 
interface fT61 . Better estimates of the effective properties can probably be obtained by incorporating damage in the 
models of PBX 9501 discussed above. However, this issue is not explored in this work because the amount and type 
of damage is difficult to quantify in polymer bonded explosives. 



5 Effective moduli of PBX 9501 

Of the various microstructures of PBX 9501 explored in the previous section, the best estimate of the Young's 
modulus of PBX 9501 is provided by the 100 particle model of the dry blend of PBX 9501 shown in Figure|8] 
overlaid by a 350x350 grid. The effective properties of PBX 9501 shown in this section have been calculated at 
different temperatures and strain rates using the 100 particle model of Figure|8] Instead of the experimentally 
determined values of moduli of HMX, data from molecular dynamics (MD) simulations 1 18 1 have been used because 
estimates using MD data of the Young's modulus of PBX 9501 at room temperature and low strain rate are slightly 
closer to the experimental value than estimates calculated using the experimentally determined modulus of HMX. 
The Young's modulus of HMX from MD simulations is 17,700 MPa and the Poisson's ratio is 0.21 compared to a 
Young's modulus of 15,300 MPa and a Poisson's ratio of 0.32 from experimental data. The Young's modulus of 
HMX is assumed constant at the temperatures and strain rates considered. 

The PBX 9501 binder data collected from various sources [19 20 6 21 1 are listed in Table0 The data are sorted in 
order of decreasing temperature followed by increasing strain rate. As shown in Table0 the Young's modulus of the 
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binder generally increases with increase in strain rate and with decrease in temperature though some of the available 
data do not follow this trend. Since values for the Poisson's ratio of the binder are not available, a value of 0.49 has 
been used for all cases. However, the value is likely to be lower near the glass transition temperature of estane 
(around -33°C). As before, a dirty binder has been used for the grid cells not fully occupied by particles. The 
differential effective medium approximation 1 2 1 has been used to calculate the effective properties of the dirty binder 
using the moduli of HMX and the binder TableQalso shows the dirty binder properties used in the simulations. It is 
interesting to note that the Poisson's ratio of the dirty binder does not change significantly even though the Young's 
modulus can vary considerably depending on the modulus of the PBX 9501 binder 

Table|8]shows the effective properties of the model RVE computed using finite elements along with the available 
experimental data on PBX 9501. The experimental data are within 15% of the finite element estimates for 
temperatures close to 25°C and at low strain rates. Good agreement is also observed for the data near 16°C and at 
high strain rates. However, at 0°C and a strain rate of 2200/s the experimental value is around 2 times the predicted 
values of Young's modulus while at -20°C and high strain rate, the experimental value is half the estimate. 
Experimental data on the Poisson's ratio of PBX 9501 are available only at 25°C and 0.005/s strain rate. The 
predicted Poisson's ratio under these conditions is considerably lower that the experimental value. Since the 
Poisson's ratio is highly sensitive to numerical error 1221 . more accurate numerical calculations may be required to 
resolve this issue. 

These results show that, based on existing data on PBX 9501, it is not possible to conclude that the effective moduli 
of PBX 9501 can be predicted accurately if the moduli of its constituents are known for various temperatures and 
strain rates. However, the results do show that for certain RVEs, close approximations to the effective Young's 
modulus of PBX 9501 can be obtained from two-dimensional finite element simulations. 



6 Summary and conclusions 

A two-dimensional finite element based method for approximating the effective elastic moduli of particulate 
composites has been presented. This approach has been applied to representative volume elements (RVEs) containing 
circular particles for volume fractions from 0. 1 to 0.92. The estimated Young's moduli and Poisson's ratios have been 
found to closely approximate those from differential effective medium approximations. Comparisons with 
three-dimensional finite element simulations have also shown that the two-dimensional approach provides good 
approximations for RVEs containing polydisperse spherical particles. The estimates of Young's modulus from 
three-dimensional models have been found to be slightly higher than from two-dimensional models. On the other 
hand, the Poisson's ratios from three-dimensional models have been found to be lower than estimates from 
two-dimensional models. 

Since polymer bonded explosives such as PBX 9501 contain two widely separated particle sizes, a number of 
manually generated microstructures of PBX 9501 have been simulated. The models were designed to contain a few 
large particles and a larger number of smaller interstitial particles while avoiding particle-particle contact. Such 
models have been found to underestimate the Young's modulus of PBX 9501 considerably when discretized with 
triangular elements. However, when the same models are discretized with square elements, considerably different 
estimates of Young's modulus are obtained depending on the degree of discretization. Hence, it is extremely difficult 
to determine an appropriate RVE for high modulus contrast and high volume fraction particulate composites that has 
the optimal size, number of particles, particle shapes and size distributions. The RVE also has to be such that it can be 
easily discretized using elements that are accurate and relatively computationally inexpensive. 
Randomly generated microstructures based on the particle size distribution of PBX 9501 show that the size of the 
RVE plays a lesser role in the determination of the effective Young's modulus than the amount of discretization used, 
even for highly refined meshes. If aligned square particles are used instead of circular particles, the predicted moduli 
are considerably higher than the experimental moduli because contacts between particles increase dramatically. It is 
also observed that the numerical simulations consistently underestimate the effective Poisson's ratio. Since the 
Poisson's ratio is highly sensitive to numerical error, accurate numerical simulations have to be performed to resolve 
this issue. 

Simulations of PBX 9501 for various temperatures and strain rates suggest that if a RVE is chosen that predicts a 
reasonable value of Young's modulus at a given temperature and strain rate, it can be used to obtain acceptable 
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estimates at other temperatures and strain rates. Further experimental data are required to confirm this finding. 
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A Two- and Three-Dimensional Moduli 

The determination of the effective properties of random particulate composites using numerical techniques requires 
the use of a three-dimensional representative volume element (RVE). If a two-dimensional RVE is selected, any plane 
strain representation immediately implies that the material is composed of unidirectional cylindrical particles. On the 
other hand, if a plane stress representation is used, the physical interpretation of the composite is a thin sheet cut from 
the particulate composite. 

One of the intents of this work has been to show that acceptable estimates of the effective elastic moduli of particulate 
composites can be obtained from two-dimensional plane strain calculations of the stress and strain fields. However, 
two-dimensional calculations in a single plane of the composite lead to two-dimensional or planar effective elastic 
moduli. For example, if a uniform uniaxial displacement ui is applied to a unit square of an isotropic, homogeneous, 
linear elastic material with Young's modulus E and Poisson's ratio v in plane strain, the resulting reaction force is Fi 
and the displacement perpendicular to the direction of the applied displacement is —U2- Since the square is of unit 
size, the nonzero strains are en = ui and £22 — ~-U2- The nonzero stresses are an — Fi and 0-33 = vFi, the latter 
obtained from the constitutive relation 

= {l/E)[{l + v)ai.j - vaj,k5,j] i, j,k ^ 1,2,3. (6) 

Note that E and i' are three-dimensional or true linear elastic moduli. It can be shown that the apparent Poisson's 
ratio in the plane (e22/eii) is given by 

i^2D = (22/^11 = U2/U1 = v/{l - ly) (7) 

and the apparent Young's modulus in the plane (crii/eii) is given by 

E2D = fTii/eii = F^/ui = E/{1 - 1.2) (8) 

These plane moduli, E2D and j/2d, are referred to as two-dimensional moduli following Jun and Jasiuk ^SJ and 
Torquato 1231 . The reason for such a name is that these apparent moduli are the same as those obtained if the 
constitutive equations are expressed, in two-dimensional form, as 

^ {1/E2d)[{1 + i^2d)o-ij - iy2D<ykkSij] i,j,k^ 1,2. (9) 

The results presented in this work suggest that these relations between two- and three-dimensional elastic moduli 
(which are exact only for isotropic, homogeneous, linear elastic materials) can be used to obtain reasonable 
approximations of three-dimensional elastic moduli of random particulate composites from two-dimensional 
numerical simulations. 
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Figure captions 



Figure 1. RVEs containing 10% to 92% by volume of circular particles, fp is the volume fraction of particles in a 
RVE. 

Figure 2. Comparison of finite element (FEM) and differential effective medium (DEM) predictions. Eb is the 
Young's modulus of the binder. 

Figure 3. Two- and three-dimensional RVEs containing 70%, 75% and 80% particles by volume, fp is the volume 
fraction of particles. 

Figure 4. Two- and three-dimensional estimates of effective Young's modulus and Poisson's ratio. 

Figure 5. Microstructure of PBX 9501 (adapted from fTl). The largest particles are around 300 microns in size. 

Figure 6. Manually generated microstructures containing ^ 90% by volume of circular particles. 

Figure 7. Microstructure containing 92% particles modeled with six-noded triangles and four-noded squares. 

Figure 8. Circular particle models based on the dry blend of PBX 9501 . 

Figure 9. Circular particle models based on the particle size distribution of pressed PBX 9501. 

Figure 10. Square particle model based on the size distribution of pressed PBX 9501. 
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Tables 
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Table 1: Effective moduli of the six model PBX 9501 microstructures from FEM calculations six-noded triangular 
elements. E is the Young's modulus and v is the Poisson's ratio. 





Expt. 






Model RVE 






Mean 


Std. Dev. 






1 


2 


3 4 


5 


6 






E (MPa) 


1013 


116 


126 


130 42 


183 


192 


132 


54 


v 


0.35 


0.34 


0.32 


0.32 0.44 


0.28 


0.25 


0.33 


0.07 



15 



Table 2: Effective elastic moduli of the four models of the dry blend of PBX 9501 . 



No. of Size Young's Modulus (MPa) Poisson's Ratio 

Particles (mm) FEM Expt. FEM Expt. 

256x256 350x350 256x256 350x350 

loo 065 1959 968 ToB 022 O20 035" 

200 0.94 2316 1488 1013 0.23 0.23 0.35 

300 1.13 2899 2004 1013 0.25 0.24 0.35 

400 1.33 4350 2845 1013 0.25 0.25 0.35 
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Table 3: Material properties of the dirty binder for the four pressed PBX 9501 microstructures. fp is the volume 
fraction of particles. 



No. of 


fp 




Properties of dirty binder 


Particles 




in binder 


E (MPa) 


i/ 


100 


0.89 


0.28 


1.583 


0.484 


200 


0.87 


0.37 


2.1395 


0.481 


500 


0.86 


0.43 


2.713 


0.478 


1000 


0.855 


0.45 


2.952 


0.477 
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Table 4: Effective elastic moduli of the four models of pressed PBX 9501. 



No. of Size Young's Modulus (MPa) Poisson's Ratio 

Particles (mm) FEM Expt. FEM Expt. 

256x256 350x350 256x256 350x350 

loo 036 2925 1798 ToB 023 OAA 035 

200 0.42 3342 2408 1013 0.21 0.20 0.35 

500 0.54 5256 3994 1013 0.25 0.23 0.35 

1000 0.68 6171 4756 1013 0.25 0.25 0.35 
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Table 5: Material properties of the dirty binder for the three PBX microstructures with square particles. 



No. of 


fp 




Properties of dirty binder 


Particles 




in binder 


E (MPa) 


i> 


700 


0.868 


0.39 


2.358 


0.4799 


2800 


0.866 


0.40 


2.448 


0.4795 


11600 


0.863 


0.42 


2.588 


0.4788 
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Table 6: Effective stiffnesses of the model microstructures with square particles. 



No. of 


Size 


Young's Modulus (MPa) 


Poisson's Ratio 


Particles 


(mm) 


FEM (256x256) Expt. 


FEM (256x256) Expt. 


700 


3.6 


9119 1013 


0.26 0.35 


2800 


5.3 


9071 1013 


0.27 0.35 


11600 


9.0 


9593 1013 


0.27 0.35 
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Table 7: Material properties used to calculate the effective moduli of PBX 9501 for different temperatures and strain 
rates. Numbers in square brackets refer to the source of the binder data. 



Temp. 


Strain 


HMX 


PBX 9501 Binder 


'Dirty' Binder 




Rate 


Young's 


Poisson's 


Young's 


Poisson's 


Young's 


Poisson's 






Modulus 


Ratio 


Modulus 


Ratio 


Modulus 


Ratio 


(°C) 


(/s) 


(MPa) 




(MPa) 




(MPa) 




25 


0.005 


17700 


0.21 


0.59 L6J 


0.49 


2.15 


0.479 


25 


0.008 


17700 


0.21 


0.73 m 


0.49 


2.66 


0.479 


25 


0.034 


17700 


0.21 


0.81 


0.49 


2.95 


0.479 


25 


0.049 


17700 


0.21 


0.82 


0.49 


2.99 


0.479 


25 


2400 


17700 


0.21 


300 EOl 


0.49 


1001 


0.474 


22 


0.001 


17700 


0.21 


0.47 ED 


0.49 


1.71 


0.479 


22 


1 


17700 


0.21 


1.4 ED 


0.49 


5.10 


0.479 


22 


2200 


17700 


0.21 


3.3 ED 


0.49 


12.0 


0.479 


16 


1700 


17700 


0.21 


22.5 fT9l 


0.49 


81.5 


0.479 





0.001 


17700 


0.21 


0.85 I21J 


0.49 


3.10 


0.479 





1700 


17700 


0.21 


2461191 


0.49 


833 


0.475 





2200 


17700 


0.21 


4ED 


0.49 


14.6 


0.479 


-15 


0.001 


17700 


0.21 


1.4 ED 


0.49 


5.10 


0.479 


-15 


1 


17700 


0.21 


5.7 ED 


0.49 


20.7 


0.479 


-15 


1000 


17700 


0.21 


16001^ 


0.49 


4082 


0.458 


-20 


0.001 


17700 


0.21 


1-6 [2LI 


0.49 


5.83 


0.479 


-20 


1200 


17700 


0.21 


1600 ED 


0.49 


4082 


0.458 


-20 


1700 


17700 


0.21 


1333 1191 


0.49 


3556 


0.461 


-40 


0.001 


17700 


0.21 


5.7 LZU 


0.49 


20.7 


0.479 


-40 


0.001 


17700 


0.21 


5.3 lT9l 


0.49 


19.3 


0.479 


-40 


1300 


17700 


0.21 


1000 f2T1 


0.49 


2838 


0.464 
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Table 8: Effective Young's modulus from finite element simulations of a 100 particle model of the dry blend of PBX 
9501. Numbers in square brackets refer to the source of the PBX 9501 data. 



Temp. 


Strain 


FEM Estimates 


Experiments 




IVtlLC 






vV^ n M rr c 








/^/nn 1 n c 


IvaLlU 




Iva LIU 






V^iVll a. ) 




a. ) 










n 1 Q 


1 (\A(\ 1 9/1 1 ^ 






U.Ui 




yj.zXf 


11C\ \1C\\ 
1 l\J \zXj \ 

1020 1191^ 




25 


0.034 


1123 


0.20 






25 


0.05 


1129 


0.20 


1013 


0.35 


25 


2400 


11273 


0.31 






22 


0.001 


926 


0.17 






22 


1 


1399 


0.23 






22 


2200 


2050 


0.26 






16 


1700 


4981 


0.30 


4650 fT9l'= 







0.001 


1144 


0.21 









1700 


10808 


0.31 









2200 


2243 


0.27 


54801191 




-15 


0.001 


1399 


0.23 






-15 


1 


2652 


0.28 






-15 


1000 


14602 


0.28 






-20 


0.001 


1481 


0.23 






-20 


1200 


14602 


0.28 






-20 


1700 


14291 


0.29 


6670 1113'^ 




-40 


0.001 


2652 


0.28 






-40 


0.001 


2562 


0.27 






-40 


1300 


13778 


0.29 


12900 ^191'' 





a - experimental data are for a strain rate of 0.001/s. 

b - experimental data are for a strain rate of 0.01 1/s and a temperature of 27°C. 
c - experimental data are for a strain rate of 2250/s and a temperature of 17°C. 
d - experimental data are for a strain rate of 2250/s. 
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Figures 
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Particle Vol. Frac. Particle Vol. Frac. 

(a) Young's Modulus. (b)Poisson's Ratio. 

Figure 2: Comparison of finite element (FEM) and differential effective medium (DEM) predictions. is the Young's 
modulus of the binder. 
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fp - 0.7 fp = 0.75 U = 0-8 




(a) Two-dimensional RVEs. 
fp = 0.7 fp = 0.75 fp = 0.8 




(b) Three-dimensional RVES. 

Figure 3: Two- and three-dimensional RVEs containing 70%, 75% and 80% particles by volume, fp is the volume 
fraction of particles. 
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Particle Vol. Frac. Particle Vol. Frac. 

(a) Young's Modulus. (b) Poisson's Ratio. 

Figure 4: Two- and three-dimensional estimates of effective Young's modulus and Poisson's ratio. 
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Figure 5: Microstructure of PBX 9501 (adapted from |7J). The largest particles are around 300 microns in size. 
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Model 1 Model 2 Model 3 




Model 4 Model 5 Model 6 

Figure 6: Manually generated microstructures containing ~ 90% by volume of circular particles. 
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(a) 92% particles. 


(b) Model 6 - 90% particles. 






^^^^^^^^^^^^^^^^ i 

(c) Six-noded Triangles. 


(d) Four-noded Squares 


■ 





Figure 7: Microstructure containing 92% particles modeled with six-noded triangles and four-noded squares. 
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0.65x0.65 mm2 0.94x0.94 mm^ l.Bxl.lSmm^ 1.33 x 1.33 mm^ 

100 Particles 200 Particles 300 Particles 400 Particles 

Figure 8: Circular particle models based on the dry blend of PBX 9501 . 
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0.36x0.36mm2 0.42x0.42 mm^ 0.54x0.54 inm^ 0.68x0.68 mm^ 

100 Particles 200 Particles 500 Particles 1000 Particles 

Figure 9: Circular particle models based on the particle size distribution of pressed PBX 9501. 
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3.6x3.6 mm^ 5.3x5.3 inm^ Q.OxQ.Omm^ 

700 Particles 2800 Particles 1 1600 Particles 

Figure 10: Square particle model based on the size distiibution of pressed PBX 9501. 



33 



